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ABSTRACT 

We examine the long-standing cooling flow problem in galaxy clusters with 3D MHD simulations of 
isolated clusters including radiative cooling and anisotropic thermal conduction along magnetic field 
lines. The central regions of the intracluster medium (ICM) can have cooling timescales of ~ 200 Myr 
or shorter — in order to prevent a cooling catastrophe the ICM must be heated by some mechanism such 
as AGN feedback or thermal conduction from the thermal reservoir at large radii. The cores of galaxy 
clusters are linearly unstable to the heat-flux-driven buoyancy instability (HBI), which significantly 
changes the thermodynamics of the cluster core. The HBI is a convective, buoyancy-driven instability 
that rearranges the magnetic field to be preferentially perpendicular to the temperature gradient. For 
a wide range of parameters, our simulations demonstrate that in the presence of the HBI, the effective 
radial thermal conductivity is reduced to < 10% of the full Spitzcr conductivity. With this suppression 
of conductive heating, the cooling catastrophe occurs on a timescale comparable to the central cooling 
time of the cluster. Thermal conduction alone is thus unlikely to stabilize clusters with low central 
entropies and short central cooling timescales. High central entropy clusters have sufficiently long 
cooling times that conduction can help stave off the cooling catastrophe for cosmologically interesting 
timescales. 

Subject headings: convection — galaxies: clusters: general — instabilities MHD — plasmas — X-rays: 
galaxies: clusters 



1. INTRODUCTION 

X-ray observations of the intracluster medium (ICM) 
of relaxed galaxy clusters show a centrally peaked sur- 
face brightness distribution. The observed temperatures 
and densities are high enough that the p lasma can hav e 
a cooling time much less than 500 Myr (jSarazi n 1986). 
The standard isobaric cooling flow model predicts mass 
dropping out of the X-ray emitting ICM at rates in ex- 
cess of 100-500 M yr _1 . However, X-ray spectroscopic 
observations with Chandra and XMM-Newton have ruled 
out classical cooling flows of material below 1 keV as it 
would co piously emit lines such as Fe XVII that are not 
observed (|Peterson fe Fabianll2006f h Therefore, a mech- 
anism is required to heat the ICM to avert this cooling 
catastrophe. 

The plasma in the ICM has temperatures of 1-15 keV 
and number densities of 10 _3 -10 _1 cm' 3 . The mag- 
netic field in the ICM is estimated to be in the range 
of 0.1 -10 nG depending on where the measurement is 
made (jCarilli fe Ta ylor 2002). Under these conditions, 
the Coulomb mean free path is many orders of magni- 
tude larger than the gyroradius; e.g., at T = 3 keV, 
n e = 10 _2 cm -3 , and B = 1 /iG, the mean free path is 
A m fp ~ 0.3 kpc, while the electron gyroradius is p e ~ 10 s 
cm. The mean free path is, however, smaller than the 
temperature gradient scale length. As a result, a fluid de- 
scription of the ICM plasma, e.g. MHD, is appropriate, 
but the effects of anisotropic heat and momentum trans- 
port must also be incl uded. The Braginskii-MHD equa- 
tions ( Bragins kiill965h are a modification of the standard 
MHD equations to include anisotropic transport due to 
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the magnetic field. 

AGN feedback and conduction from the thermal bath 
at large radii are two of the most often discussed mech- 
anisms for heating cool cluster cores. This work shall 
only briefly consider the former. The latter, ther- 
mal conduction, has been studied by many authors 
("e.g.. iBinnev fe Cowidl981tlNaravan fe Medvedev>20Qll : 
IZakamska fe Naravanll2003HGuo et al.ll2008fh Because of 
uncertainties associated with the magnetic geometry, the 
heat flux is often parameterized as an effective thermal 
conductivity given as a fraction, /g p , of the ideal (field- 
free) Spitzer heat flux. Previous work, however, has not 
considered the dynamic consequences of anisotropic con- 
duction. The plasma in galaxy clusters is unstable to two 
different convective instabilities driven by anisotropic 
thermal conduction along magnetic field li nes. The first , 
the magnetotherm al instability, or MTI (jBalbusI 120001 : 
IParrish et all [20081. occurs when the temperature gra- 
dient and gravity are in the same direction, as is true 
at large radii in galaxy clusters. Well inside the cooling 
core, the heat-flux-driven-buoyancy instability, or HBI, 
occurs where the te mperature gradient is in the opposite 
direc tion of gravity (|Quataertj|2008l : IParrish fe Quataertl 
120081) . The HBI occurs in the cooling core of the ICM 
where the temperature increases outward. Nonlinear 
simulations of the MTI and HBI have shown that they 
significantly modify the thermal conductivity, because 
they saturate by rearranging the magnetic field geometry. 
Thus, determining the effective conductivity in the ICM 
requires considering the back-reaction of the anisotropic 
heat flux on the magnetic field geometry. 

In this work, we examine the stability of the cooling 
cores of galaxy clusters using three-dimensional MHD 
simulations including anisotropic thermal conduction 
and cooling. In particular we assess the interplay be- 
tween cooling and the HBI. This paper is organized as 
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follows. In ij2]we summarize the physics of the HBI and 
the thermal instability and their nonlinear saturation in 
local calculations. In fJ3]we then describe the equations 
of Braginskii-MHD and the numerical tools we utilize to 
solve them. §4] presents our fiducial cluster model (based 
on Abell 2199) in detail and examines its evolution over 
cosmic time. In fJSJ we examine a variety of variations 
of our fiducial model including cluster halo parameters, 
magnetic field strength and geometry, and the central 
cluster entropy. We also show a few experiments with 
a simple AGN heating model ( £|5.7p . Finally, in S|6] we 
discuss the implications of this work for the cooling flow 
problem and highlight some of our plans for future work. 

In parallel t o ou r work described here, 
iBogdanovic et all (|2009f ) has conducted a similar 
study using similar numerical methods. Their simula- 
tions start with different initial conditions and cover 
a different part of parameter space but reach broadly 
similar conclusions. 

2. PHYSICS OF THE HBI AND COOLING 

The linear physics of the HBI and its n onlin- 
ear evolution are outlin ed in iQuataert] (2008) and 
iParrish fe Quataertl (|2008| ). respectively. We briefly re- 
view them here for clarity. The heat-flux-driven buoy- 
ancy instability is a convective instability driven by a 
background heat flux with the temperature gradient as 
the source of free energy. In contrast, the entropy gra- 
dient drives the more familiar Schwarzschild convection. 
For an arbitrarily orie nted magnetic fi eld, the dispersion 
relation of the HBI is (jQuataertll2008h : 



0=LOU 2 



- 2 

" IWcond^ 



dlnT 

dz 



k 2 

K) k 2 



(1 



where N is the Brunt- Vaisala frequency, 

1 dPdlnS 



N = -- 



and 



oj 2 



7P dz dz 

uj 2 ~(k-v A ) 2 , (3) 

where va = B / (inp) 1 / 2 is the Alfven speed, and 

2 /- \ 2 
Wcond = [b ■ k) (4) 

is the frequency for conduction to act on a given scale, 
with b the unit vector directed along the magnetic field 
and x the thermal diffusivity 3 in units of cm 2 s _1 . This 
dispersion relation is written without loss of generality 
for a geometry in which gravity and the initial atmo- 
spheric gradients are in the ^-direction and the initial 
magnetic field lies in the x-z plane. 

For a weak magnetic field, equation ([I]) has unstable 
solutions for either sign of the temperature gradient. The 
case of dT/dz > corresponds to the HBI. For a weak, 
initially vertical magnetic field (b z = 1, b x = 0), the 
growth rate of the HBI is given by 



dlnT\ k\ 



dz 



k 2 ' 



(5) 



3 The literature is not consistent regarding the use of x and «• 
We will use \ to represent a true diffusion coefficient and k to 
represent a conductivity in units erg cm -1 s _1 K — . 



where k± is with respect to gravity. Qualitatively, one 
can picture the HBI as being driven by regions of con- 
verging and diverging perturbed magnetic field lines. 
Regions of converging magnetic field are conductively 
heated and become buoyant. In local simulations, the 
HBI generates MHD turbulence and a modest magnetic 
dynamo that amplifies the field. The most prominent 
method by which the HBI saturates is via a significant 
reorientation of the magnetic field geometry. The HBI 
takes a largely vertical field and reorients it to become 
largely horizontal. This fact is crucial for cluster cores 
since this reorientation of the magnetic field can signifi- 
cantly reduce the heat transport across an HBI unstable 
region. 

In addition to driving the HBI, thermal conduction 
also has a significa nt impact on th e rmal instability as 
originally shown in iField fc Saslawl (|1965l ). The Field 
criterion states that wavelengths longer than 



A, 



Tk 



n e n p A(T) 



1/2 



(6) 



are thermally unstable, where A(T) is the cooling func- 
tion discussed later. For modes with wavelengths shorter 
than Xp, the conduction time is shorter than the cool- 
ing time, and local perturbations are stabilized. In the 
ICM, the plasma is often locally stable to thermal in- 
stability, but unsta ble to global modes ev en in the pres- 
ence of conduction (Ki m fc Naravanl feOOS) . Equation © 
was derived under the assumption of isotropic conduc- 
tion. The results are similar for anisotropic conduction, 
except that conduction can only stabilize perturbations 
with short wavelengths along the magnetic field. 

An interesting way to examine the p hysics of cooling i n 
clusters comes from the recent work of lVoit et al.l ((2008) , 
W who examined the role of the central entropy of the ICM 
as an indicator of feedback and star formation in galaxy 

clusters. The en tropy is defined as K = ksTrie 2 ^ 3 . In 
(2) iVoit et all (|2008f ). the entropy profile for a cluster is fit 



using 



K(r) = K a + K 



too 



100 kpc 



(7) 



where Kq is approximately the central entropy, .Kioo is 
the power law normalization and a > is the power 
law exponent. They find that low entropy clusters, those 
with Kq < 30 keV cm 2 , have stronger Ha emission, star 
formation indicators, and AGN activity than higher en - 
tropy clusters, K > 30 keV cm 2 l|Cavagnolo et al.ll2008f) . 
As a matter of terminology, clusters with an inwardly 
decreasing temperature are referred to as cool-core or re- 
laxed clusters. Clusters with an isothermal or inwardly 
increasing temperature profile are referred to as non cool- 
core clusters. 

These observational results can be qualitatively under- 
stood in light of the Field criterion (eq. [6]). When cool- 
ing is pure Bremsstrahlung, the Field length becomes a 

3/2 1/2 

function of entropy, scaling as Af oc K f s ' . Thus, 
cooling can take place on small length scales when the 
entropy is low. The HBI decreases /s p , making smaller 
wavelengths unstable to cooling. Fully non-linear simu- 
lations, such as those presented here, are needed to un- 
derstand the combined dynamics of cooling and the HBI. 
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3. METHOD 

3.1. Equations of MHD with Anisotropic Heat 
Conduction 

We solve the usual equations of ideal MHD with the 
addition of a vector heat flux, Q, and a gravitational 
acceleration g = — V$: 



! + v.( H = o, 



d(pv) 
dt 





( B 2 > 




pvv + 




i'-71 



dE 
~dt 



-V- 
V Q 

dB 

~dt 
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8tt 



B{B-v) 



-in 



h pV$ -v = H-£, 
- V x (v x B) = 0, 



(8) 

= 0, 
(9) 

(10) 
(11) 



where the symbols have their usual meaning. The total 
energy E is given by 



E 



B B 



2 8tt 
= p/( 7 -l) 



and the internal energy, e 
5/3 throughout. 
The anisotropic heat flux is given by 



(12) 

We assume 7 = 



Q = -nksXc(T, n)bb • VT, 



(13) 



where the thermal diffusivity is given by the Spitzer value 
(|Spitzerl 11962]) and b is a unit vector in the direction of 
the magnetic field. The Spitzer thermal diffusivity is 
given by 



X c(T,n) = 8xl0 31 



T \ 5/2 

lOkevJ 



5 x IO- 3 



9 — 1 

cm s . 



(14) 

Note that % is a thermal diffusivity, and it depends 
inversely on the density. The Spitzer conductivity is 
K Sp = n ^BXc which has only the well-known T 5 / 2 de- 
pendence and no density dependence. 

The energy equation also includes heating (Ti.) and 
cooling (£) terms . The cooling function we adopt is from 
iTozzi fc Norman] (|2001f ) with the functional form 



n e n p A(T), 



(15) 



with units of erg cm~ 3 s . The temperature dependence 
is a fit to cooling dominated by Bremsstrahlung above 1 
keV and metal lines below 1 keV with 



A(T) = [diksT)- 1 - 7 + C 2 (k B T) - 5 + C 3 ] 10 



-22 



(16) 



where d = 8.6 x 10~ 3 , C 2 = 5.8 x 10~ 2 , and C 3 = 
6.3 x 10~ 2 , for a metallicity of Z — 0.3Z©, with units 
of [Ci] = ergcm 3 s _1 . The majority of our runs do not 
include additional heating terms. For these runs H = 
in equation (flQ)) . 

In runs with heating, we adopt a heating profile of the 
form 

■H = H e- (r/rH) \ (17) 



TABLE 1 

TlMESCALES IN MODEL Tl 



Timescale 


Symbol 


Time a (Myr) 


Sound Crossing 


TS 


45 


Alfven Crossing (1/iG) 


TA 


1.1 x 10 3 


Conduction 


^cond 


20 


HBI Growth 


thbi 


120 


Cooling 


''"cool 


1.4 x 10 3 



a Evaluated for L = 50 kpc as volume-averaged 
quantities. 

where rn is the scale radius of heating and the normal- 
ization is chosen as 



Ho 



(18) 



H' 



where Lthcrm is the total thermal heating input. This 
model is motivated by a simple description of AGN feed- 
back. We discuss simulations with heating in more detail 
in fT3 

3.2. Timescales 

We now discuss several of the key timescales in galaxy 
clusters in order to provide some intuition for the impor- 
tant physical processes (See Tabled]). We examine these 
timescales for volume-averaged quantities in our fiducial 
atmosphere of A2199. The generation of the fiducial 
atmosphere is discussed in §4.11 The volume-averaged 
sound speed is approximately 10 8 cm s , corresponding 
to a sound crossing time over 50 kpc of t$ ~ 45 Myr. 
A weak magnetic field of 1 nG corresponds to an Alfven 
crossing time over 50 kpc of ta ~ 10 12 years. 

For a more typical magnetic field of 1 /iG, the Alfven 
speed is 4.4 x 10 6 cm s _1 , leading to a crossing time across 
50 kpc of 1.1 Gyr. We will see shortly that the Alfven 
timescale is the longest timescale in the problem. For this 
cluster the central magnetic beta, the ratio of thermal 
pressure to magnetic pressure is /3q = 8np/B 2 « 6,600 
for B = 1 /iG. Even for B = 10 fiG, the central magnetic 
beta is significantly greater than unity. Further out in 
the core, where the thermal pressure has dropped, the 
beta parameter is typically of order several hundred. 

Also of interest is the heat conduction timescale. Our 
model cluster has a volume-averaged thermal diffusiv- 
ity of (x) ~ 3.8 x 10 31 cm 2 s _1 , which yields the scale- 
dependent conduction time 

L 2 



X 



20 Myr (L 
0.80 Myr {L -- 



50 kpc) 
10 kpc). 



(19) 



The HBI growth time in the fast conduction limit is 



THBI 



dlnTd0 
dr dr 



-1/2 



126 Myr. 



(20) 



As we will see later, the HBI has an opportunity to grow 
significantly compared to the average time between ma- 
jor mergers, roughly 5 Gyr, for a typical cluster. The 
final timescale of interest is the cooling time at the cen- 
ter of the cluster 

7 e 



Tcool — 



1.4 Gyr. 



(21) 



7-1 n e n p A(T) 

This cooling time estimate is in general too long (by al- 
most a factor of 2) , as it does not account for the increase 
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in the cooling rate as temperature decreases and the den- 
sity and line emission increase. Nonetheless, the cooling 
time is much longer than the HBI growth rate for the 
fastest growing, short-wavelength modes. Thus, the HBI 
growth is likely to play a significant role in the thermal 
evolution of the cluster core. 

3.3. Numerical Tools 

For our simulations we use the Athena MHD 
code dGardiner fc Stone! 120081 : IStone et alj l2008fl com- 
bined with the anis otrop ic conduction methods o f 
iParrish fc Stone! (|2005t) and lSharma k. Hammettl (|2007l ). 
In particular, we use harmonic averaging of the conduc- 
tivity and the monotonized central difference limiter on 
transverse heat fluxes to ensure stability. The heating, 
cooling, and thermal conduction are operator split and 
sub-cycled with respect to the MHD timestep. The cool- 
ing simulations are implemented with a temperature floor 
of T = 0.05 keV, below which UV lines become impor- 
tant, and the cooling curve fit is no longer accurate. This 
temperature floor prevents the cooling catastrophe from 
going to completion. 

Most of the simulations in this work are done on uni- 
form Cartesian grids of (128) 3 . One high-resolution run 
is done at (256) 3 . We use modified reflecting boundary 
conditions for all MHD variables, in which the pressure 
and density are extrapolated in the ghost zones, but the 
other variables are reflected. This prevents the gravita- 
tional source term from introducing a spurious acceler- 
ation at the boundary. For the temperature boundary 
condition, we introduce a thermal bath at r > r max with 
a fixed temperature T outcr - For almost all of our runs 
fmax = 200 kpc. This thermal bath physically repre- 
sents the massive thermal energy available in the ICM 
outside the cluster core. Since we are simulating the en- 
tire cluster in a Cartesian geometry, there is no inner 
boundary condition at the cluster center. 

To seed multiple modes of the HBI and break symme- 
try, we add Gaussian white noise perturbations to the 
initial velocity field such that the applied perturbation is 
everywhere a fixed fraction of the sound speed, typically 
w 1%. In the absence of the HBI, cooling, or an ap- 
plied perturbation, we find that we can hold hydrostatic 
equilibrium to better than one part in 10 4 . 

All of our runs, unless otherwise noted, are run to a 
time of 9.5 Gyr, a large fraction of the age of the universe. 
In a small number of runs, the simulations do not go 
to completion as a result of very severe cooling flows 
concentrating large quantities of mass into the central 
few zones of the grid. These exceptions are noted. 

4. FIDUCIAL SIMULATION 

4.1. Initial Conditions 

To introduce the phenomenology of the HBI in cluster 
cores, we begin by discussing a simple fiducial calcula- 
tion based on observ ations of the galaxy c l uster Abell 
2199 as discussed in p akamsk a~fc Naravanl (|2003l ) and 
Uohnstone et al.l ((2002). This fiducial run is identified 
as Tl in the table of Runs (Table [2]). Our initial con- 
ditions for the cluster core are obtained by constructing 
a spherically symmetric atmosphere in both hydrostatic 
and thermal equilibrium. We have found that it is quite 
advantageous to start from thermal equilibrium. Runs 



without thermal equilibrium experience thermal fronts 
propagating from the boundaries. Starting in thermal 
equilibrium produces a much more physical initial con- 
dition. The equations for this equilibrium are 



dP 
dr 



knT dr 



1 d 



dr 



dr 



72 fsp K ~ J = n e n p k{T) - H, 



(22) 



(23) 



where fs p is the initial effective thermal conductivity and 
the heating function is neglected for our fiducial case. 
Our potential is chosen to be a softened NFW potential 
with a dark matter density distribution given by 



PDA! 



M /2tt 



(r + r c )(r + r s ) 2 ' 



(24) 



where Mq is the scale mass, r s is the scale r adius, and r c 
is the softening radius (|Navarro et al.lll997l ). The poten- 
tial softening is important for numerically maintaining a 
very accurate hydrostatic equilibrium. The potential is 
given by 



>=-2GM 



In 



1 + r/r c 
1 + r/r s 



In (1 + r/r c ) 



r/r c 



r s (r s - 2r c ) ln(l + r/r s ) 



r c (r s ~ r c ) 2 



r /r c 



(25) 



The plasma is modeled as a fully ionized ideal gas with 
(i = 0.617 and fi e = 1-176. 

We solve equations (|22|) and (|23|) as a two point bound- 
ary value problem with the constraints of matching T at 
both the inner (Tj) and outer (T ) boundary. As these 
are a third-order system of equations, we choose a further 
symmetry constraint, namely, that the heat flux vanishes 
at the center. This system of equations is slightly stiff, 
but generally is soluble with a shooting method with a 
good initial guess. We only solve these ODEs to estab- 
lish our initial condition. The subsequent evolution of the 
equilibrium is calculated using the full system of partial 
differential equations (PDEs) for MHD, equations ((HJ- 

(ED- 

For our fiducial model, we use a physical initial mag- 
netic field geometry, that of tangled magnetic fields. 
First, in constructing the two-point boundary value for 
our initial conditions, we set /s p = 1/3. Second, the 
initial fields are tangled with a Kolmogo rov spectrum 
using the method outlined in detail §4.2 of IParrish et al.l 
(2008). We initialize the fields in Fourier space as 



A(k) = A 



&peak 



(26) 



where fc pea k is chosen as the wavenumber corresponding 
to 2-4 times the grid scale. We also choose a low-A; cut- 
off corresponding to wavelengths of 1/2 to 1/4 of the 
domain size. We randomize the phase and use the Fast 
Fourier Transform to calculate the vector potential in 
real space. We choose a = —17/6, the appropriate k- 
space scaling for the ID power spectrum of Kolmogorov 
turbulence. Our last step in initializing the magnetic 
field is to difference the vector potential to calculate a 
manifestly divergence-free initial field. The normaliza- 
tion of the magnetic field for our fiducial run is such that 
(\B\) = 10" 9 G. 
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TABLE 2 

Initial Properties of Nonlinear Runs 



Run 


M (M ) R s (kpc) 


-R c a (kpc) 


Ti (keV) 


T (keV) 


-Bo (G) 


K (keV cm 2 ) 


Rl 


3.8 x 10 14 


390 


20 


2 


5 


10~ 9 , radial 


15.5 


T1 b 


3.8 x 10 14 


390 


20 


2 


5 


10 -9 , tangled 


22.4 


Tl-HB 


3.8 x 10 14 


390 


20 


2 


5 


10~ 6 , tangled 


22.4 


Tl-256 


3.8 x 10 14 


390 


20 


2 


5 


10~ 9 , tangled 


22.4 


T2 


1.1 x 10 15 


520 


26 


1 


9.5 


10 -9 , tangled 


31.1 


T3 


3.8 x 10 14 


390 


20 


1 


6 


10~ 9 , tangled 


5.46 


El 


3.8 x 10 14 


390 


20 


3 


5 


10 -9 , tangled 


43.6 


E2 


3.8 x 10 14 


390 


20 


4 


5 


10 -9 , tangled 


83.1 


E2-NC C 


3.8 x 10 14 


390 


20 


4 


5 


10~ 9 , tangled 


83.1 


E2-HB 


3.8 x 10 14 


390 


20 


4 


5 


10 -6 , tangled 


83.1 


E3 


3.8 x 10 14 


390 


20 


4.5 


5 


10~ 9 , tangled 


122. 


E3-NC C 


3.8 x 10 14 


390 


20 


4.5 


5 


10 -9 , tangled 


122. 


Hl d 


3.8 x 10 14 


390 


20 


2 


5 


3.5 X 10" 6 , tangled 


14.0 


Il c 


3.8 x 10 14 


390 


20 


2 


5 


10 -9 , radial 


1.5 


Isol 


3.8 x 10 14 


390 


20 


4.5 


4.5 


10 -9 , tangled 


52.5 


Iso2 


3.8 x 10 14 


390 


20 


4.5 


4.5 


10 -9 , tangled 


154. 



a Softening radius of NFW halo (eq. [24]) 
b Fiducial Run 
c No conduction 



d Simulation includes additional heating (see i)5.7p . 
c Isotropic conduction only 




Fig. 1. — The temperature, density, and pressure of our fidu- 
cial initial condition are chosen to roughly correspond to those 
observed in Abell 2199. We initialize tangled magnetic fields for 
this simulation. 



Our model cluster, Abell 2199 has an inner tempera- 
ture of roughl y 2 keV and a tem perature of 5 keV near 
200 kpc (Johnston e et alJ 12002) . The gravitational po- 
tential is fit to an NFW profile with a scale radius of 
r s = 390 kpc, a softening radius of r c = 20 kpc, and a 
mass of M — 3.8 x 1O 14 M . Figure Q] shows the fiducial 
atmosphere that results from these parameters. The cen- 
tral density in our thermal equilibrium model is slightly 
lower than the observed value. 

4.2. Evolution of the Fiducial Case 



The evolution of our fiducial cluster model prominently 
illustrates the physics of the HBI and cooling in the 
galaxy cluster core. This evolution is best understood 
through a variety of diagnostics. First, in Figure [5] we 
show the time evolution of the magnetic and kinetic en- 
ergies in the cluster. There is a very weak magnetic dy- 
namo in the first 2 Gyr or so due to the HBI. The initial 
drop in magnetic energy is due to reconnection. After 
< 3 Gyr, the core loses central pressure support, inflow 
begins, and the kinetic energy increases. Correspond- 
ingly, the magnetic energy is amplified through flux- 
freezing leading to a maximum increase of A(B 2 ) m 2.5. 
We define the magnetic energy amplification at the final 
simulation time, it, as 



A(B 2 ) = 



2 ,_ (B 2 )(t f ) 



(27) 



where the angle brackets denote a volume average. The 
kinetic energy dominates the magnetic energy at all 
times. Figured] shows that the HBI is not a strong source 
of magnetic field amplification in cluster cores. 

The real hallmark of the HBI is the reorientation of 
the magnetic field geometry. Figure [3] shows the evolu- 
tion of the volume-averaged angle of the magnetic field 
with respect to the radial direction from its initial tan- 
gled state (9 = 60°) to a final saturated state of 74.6°. 
The angles given here and in Table [3] are volume av- 
erages over the entire cluster. The reorientation of the 
magnetic field is quite dramatic and takes place on just 
a few Gyr. Concomitantly, Figure [4] shows the evolution 
of the azimuthally-averaged radial-temperature profile. 
We typically bin the temperature into 5 kpc radial bins. 
As the magnetic geometry evolves to be more azimuthal, 
the thermal conduction from outside the core begins to 
shut off and the temperature starts to fall. At ~ 2.7 Gyr, 
the central temperature has hit the cooling floor of 0.05 
keV — an effective proxy for the cooling catastrophe. We 
do not remove gas from the grid after hitting the tem- 
perature floor, and thus, we never see a true "cooling 
flow." We define the time of the cooling catastrophe, t cc , 
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Fig. 2. — Time evolution of the volume-averaged magnetic and ki- 
netic energy in our fiducial run, Tl. A very weak magnetic dynamo 
is present. After 2.7 Gyr, the plasma has reached the temperature 
floor leading to a significant inflow and magnetic field increase due 
to flux freezing. 



Fig. 4. — Azimuthally-averaged radial temperature profiles in 
run Tl. The HBI effectively shuts off conduction leading to the 
cooling catastrophe that occurs around 2.7 Gyr. The temperature 
is fixed at 5 keV beyond 200 kpc 
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Fig. 3. — Evolution of the volume-averaged angle of the magnetic 
field with respect to the radial direction in run Tl. 8 = 0° is radial. 
9 = 60° corresponds to a random magnetic field. The magnetic 
field becomes significantly more azimuthal due to the HBI. 



to be when the average temperature of the inner bin has 
reached the temperature floor. 

We can explore the magnetic field amplification in 
slightly more detail. We bin the magnetic field into ra- 
dial bins, and then for each bin, j, we calculate the local 



R (kpc) 

Fig. 5. — Amplification of the magnetic energy (eq. [28]) in radial 
bins during run Tl. The HBI produces only a modest dynamo, 
which is most effective near ~ 100 kpc. The central amplification 
of the field at late times is primarily due to flux-freezing during the 
cooling catastrophe. 

amplification of the magnetic energy as 

±(B% , ^0- y (28) 

Figure [5] shows the amplification as a function of radius 
at several times. The HBI amplifies the energy most 
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efficiently in the middle of the core, just beyond 100 kpc. 
The high fields at late times in the central region are 
primarily due to flux-freezing as the density increases in 
the core. 

Finally, Figure [6] shows a p o st hoc calculation of the 
heat fluxes. See iParrish et alJ (|2008t ) §5.3 for a full dis- 
cussion of the heat flux diagnostics. We calculate the 
heat fluxes as a shell average within the radial range of 
r = 100 ± 40 kpc. We begin by defining a fiducial heat 
flux to be the radial flux through the shell if the conduc- 
tion were isotropic at the Spitzer value, namely: 

Q r = -nk BX c{T 7 n) — . (29) 
dr 

This value is the same as the heat flux with anisotropic 
conduction and purely radial magnetic field lines. We 
calculate the conductive heat flux and normalize to the 
fiducial value to calculate the Spitzer fraction, or effective 
conductivity, defined as 

/Sp = Qcond/Qr, (30) 

where Q CO nd is given by equation (flU]) . We also define a 
flux due to mass advection 

Qadv = k B {(n){T)(v r ) + {T)(6n6v r )) , (31) 

7-1 

where angle brackets denote shell averages. The final 
component of the heat flux is the convective heat flux 
given by 

Qco„v = k B ((n) (Sv r 6T) + (v r ) {SnST) (32) 
7-1 

+ (dnSTSv r )) , 

where S refers to the local deviation from the mean of 
a quantity, e.g. Sv r = v r — (v r ). The first terms of 
equations (|3"Tj) and (|32p are the dominant terms. 

The initial Spitzer fraction for tangled magnetic fields 
is /s p w 1/3 due to the average over the random field 
geometry. From run to run there is some variation in 
this initial value as there is mode power on scales larger 
than our averaging volume. As the HBI grows, the 
heat flux is reduced significantly, eventually saturating 
at /spitz = 0.07 (Fig. This dramatic reduction in 
heat flux leads to the cooling catastrophe. As the core 
cools and loses pressure support, the advective heat flux 
increases as mass is transported inwards. It is interest- 
ing to note that the convective heat flux is very small, 
especially during the HBI phase of the evolution. 

For comparison purposes, we also run our fiducial 
model with purely isotropic conductivity and radial mag- 
netic fields (run II in Table 1). The HBI is not present for 
purely isotropic conduction. Figure [7] shows the evolu- 
tion of the azimuthally-averaged temperature profile — 
the cluster reaches an almost isothermal temperature 
profile with no hint of the cooling catastrophe. At fixed 
pressure, the thermal instability can lead to either run- 
away heating or runaway cooling. This is easy to see by 
examining the form of the cooling term and conduction. 
In the Bremsstrahlung regime, cooling scales as T -3 / 2 at 
fixed pressure, while the conduction term scales as T 7 / 2 . 
For the case illustrated here, as the temperature is per- 
turbed upwards, conductive heating increases much more 
rapidly than cooling. Thus, the thermal runaway drives 
the cluster towards an almost isothermal profile. In this 
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Fig. 6. — The evolution of the components of the heat flux 
nor malized to the instantaneous fiducial (field-free) heat flux (eq. 
[29]) in a shell centered at 100 ± 40 kpc in run Tl. The heat flux 
is separated into conduction (solid line), convection (dashed line), 
and mass advection (dotted line). The final saturated conductive 
heat flux has a Spitzer fraction of 7%. 
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Fig. 7. — Evolution of the azimuthally-averaged temperature 
profile for run II with purely isotropic conduction. Thermal insta- 
bility leads to an almost isothermal state. 

simulation, there is a small amount of noise which ran- 
domizes the field a very small amount, hence (0 b) 7^ 
in Table O 

To summarize the fiducial case, we begin with a cluster 
core in hydrodynamic and thermal equilibrium. If ther- 
mal conduction were isotropic, the central temperature 
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TABLE 3 
Saturation of Nonlinear Runs 



Run 


A.(B Z ) 


max(0s) 


min/sp 


tcool (Gyr) a 


tec (Gyr)° 


Rl 


15.0 


65° 


0.13 


0.82 


4.0 


IT 


2.5 


74° 


0.07 


1.4 


2.7 


Tl-HB 


1.05 


74° 


0.06 


1 1 


2.7 


Tl-256 c 


5.9 


77° 


0.03 


1.4 


2.4 


T2 C 


1.5 


75° 


0.06 


1.1 


3.2 


T3 


2.5 


74° 


0.08 


0.20 


1.5 


El 


2.3 


75° 


0.07 


3.0 


3.5 


E2 


2.3 


76° 


0.05 


5.9 


5.7 


E2-NC 


2.2 


53° 


0.48 


5.9 


1.0 


E2-HB 


0.47 


74° 


0.08 


5.9 


7.0 


E3 


2.2 


77° 


0.05 


9.3 


>9.5 


E3-NC 


4.1 


50.6° 


0.55 


9.3 


5.9 


Hl c 


0.30 


72° 


0.09 


1.3 


>6.0 


11 


1.3 


24° 


TO 


0.82 


none 


Isol 


6.1 


67° 


0.17 


2.6 


1.7 


Iso2 


1.74 


76° 


0.06 


12.8 


>9.5 



a Initial cooling time at the innermost radii (~5 kpc). 

b Time of Cooling Cat astro phe, when the inner gridpoint has reached 

the cooling floor (see N4.2|l , 

c Run time is less than 9.5 Gyr. 



would latch onto the heating branch of the thermal insta- 
bility, continuing to rise. Instead, the HBI begins to act 
on a 100 Myr timescale, reorienting the magnetic field 
geometry, and reducing the effective thermal conductiv- 
ity from the outer part of the cluster. As the cluster 
center becomes denser and cooler, a thermal runaway 
proceeds, leading to a cooling catastrophe on a timescale 
comparable to the initial cooling time in the core of the 
cluster. 

5. VARIATION OF PARAMETERS 

In order to fully understand the behavior of the HBI 
in galaxy clusters, we now turn to an exploration of pa- 
rameter space. Table [2] lists the various runs in which we 
vary the cluster properties, magnetic field strength and 
geometry, and the central entropy of the cluster. The fol- 
lowing sections describe each of these experiments. 

Table [3] lists the saturation properties of these runs. 
The magnetic field amplification, A(B 2 ), is given as a 
volume- aver age over the cluster. The maximum of the 
magnetic field angle, max(0s), is the maximum in time of 
the volume-averaged magnetic angle. Likewise, min/g p , 
is the minimum in time of the shell-averaged heat flux. 

5.1. Radial Magnetic Fields 

To assess the importance of the initial magnetic field 
geometry, we consider a split monopole radial magnetic 
field such that B(r) = B c (r /r^)^ 2 sgn(z). This geometry 
is useful for illustrating the effects of the HBI on the 
equilibrium state. We choose a mean magnetic field of 1 
nG in order to minimize the effects of magnetic tension 
on the equilibrium state. The resulting initial condition 
has a slightly higher central pressure and density than 
our fiducial case. In general, the atmosphere is in the 
rapid conduction limit on all but the largest scales. The 
choice of radial magnetic fields gives conduction the best 
chance of thermally stabilizing the cluster. 

The evolution of the HBI in the radial field simulation 
is quite similar to that of our fiducial case. Initially the 
atmosphere is in thermal equilibrium with radial fields. 
Figure [8] shows the evolution of the volume-averaged 
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t (Myr) 

Fig. 8. — Volume-averaged angle of the magnetic field with re- 
spect to the radial direction for run Rl. = 0° corresponds to a 
radial magnetic field. The magnetic field geometry is significantly 
rearranged reaching a maximum of approximately 65° from radial. 
9(t = 0) 7^ due to discretization errors. 

magnetic field from the initially radial (9 w 0°) geome- 
try. The HBI mode grows rapidly on a timescale of ~ 100 
Myr and reorients the magnetic field to have (9) > 60° in 
just over 4 Gyr. Figure [9] shows the resultant evolution of 
the temperature profile. The atmosphere initially latches 
onto a heating branch of the thermal instability, peaking 
at a central temperature of just over 3 keV at 1.1 Gyr. In 
the absence of the HBI, it would continue to evolve to an 
almost isothermal state as in run II. Instead, however, 
the cluster undergoes a cooling catastrophe after 4 Gyr. 

The driver of this cooling catastrophe is easily seen 
from Figure [TUl which plots the heat fluxes as a function 
of time at 100 kpc. As the HBI rearranges the magnetic 
geometry, the Spitzer fraction plummets to /g p R* 0.13. 
Having reduced the contact with the thermal bath at 
large radii, cooling becomes dominant in the core and 
the central temperature starts to rapidly decrease. This 
cooling drives mass inflow to small radii giving rise to the 
large inward advective flux at late times. At all times, 
the convective heat flux is small compared to both the 
saturated conductive heat flux and the mixing length es- 
timate. Figure [Tl] shows the overall evolution of this 
cluster in a 2D slice taken at z = 0. As the magnetic 
field lines wrap in the azimuthal direction, the central 
temperature decrease is easily observed. 

The cooling catastrophe's vigor is driven by two com- 
plementary properties of the cooling curve. First, at fixed 
pressure, the cooling increases as the temperature de- 
clines at the center of the cluster. Second, below 1 keV, 
line emission from metals like iron and oxygen becomes 
increasingly important, scaling as C cx T~ 1,7 . Thus, once 
gas has cooled below 1 keV, the cooling is much harder to 
reverse. Observations show very few clusters with central 
temperatures below 1 keV. 
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Fig. 9. — Azimuthally-averaged radial temperature profiles in 
run Rl. After initially spending time on a heating branch, the HBI 
shuts off conduction leading to the cooling catastrophe around 3.6 
Gyr. 
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Fig. 10. — Time evolution of the components of the heat flux 
normalized to the instantaneous fiducial heat flux (eq. 1291 in a shell 
centered at 100 ± 40 kpc for run Rl. The heat flux is separated 
into conduction (solid line), convection (dashed line), and mass 
advection (dotted line). The final saturated conductive heat flux 
is ~ 13% of the field-free (Spitzer) value. 



5.2. Strong Magnetic Fields 

In the previous sections, we have demonstrated the 
evolution of the HBI in clusters for weak magnetic fields. 
Under these conditions, magnetic tension is not signifi- 
cant. We now consider a more realistic magnetic field of 



l^G. See ICarilli fc Tavlorl (|2002) for a review of cluster 
magnetic field measurements. When the tension force be- 
comes comparable to the buoyancy time, the HBI growth 
is suppressed at small scales. For example, for B — 1 /iG 
in our fiducial cluster, magnetic tension becomes impor- 
tant on scales smaller than 5-20 kpc depending on the 
location in the cluster. 

The cluster evolution of this case (labeled run Tl-HB) 
is similar to the weaker magnetic field tangled case. The 
maximum field strength we can simulate in this constant 
field model is limited since a modest field at the center 
can be dynamically important and low-/3 several scale 
heights out from the center. For run Tl-HB we find 
that there is a modest amount of numerical reconnection 
that dissipates some of the initial magnetic energy. A 
very weak dynamo leads to a maximum magnetic energy 
increase of only 5%. The maximum magnetic field angle 
and minimum heat flux are quite similar to the lower 
field case. In addition, the cooling catastrophe occurs at 
almost exactly the same time. Thus, a constant 1 /iG 
magnetic field provides very little stabilization for our 
fiducial cluster. 

It is interesting to examine the magnetic field geometry 
in this simulation from a different perspective, namely, 
how the average magnetic field angle varies versus ra- 
dius. This is shown in Figure [T2] Initially, the angle is 
distributed as a random variable about 60°. At 1.3 Gyr, 
the HBI has had somewhat more than 10 growth times to 
significantly increase the average angle and decrease the 
conductivity. By 2.2 Gyr, radial infall from inhomoge- 
neously cooling material has straightened the magnetic 
field out within 30 kpc, e nhancing the radial heat fl ux. 
This idea was suggested bv lBalbus &; Reynolds! (|2008| ) as 
a mechanism for opposing the HBI and slowing the cool- 
ing catastrophe. Unfortunately, by the time the magnetic 
field has been straightened out, the plasma within 20 kpc 
has a mean temperature of 0.2 keV and a mean density of 
(n e ) w 1.5 cm~ 3 . This corresponds to an increase in cool- 
ing from the initial equilibrium of C/Co ~ 350. Mean- 
while, the average angle in this region is about 50°, im- 
plying that the effective conductivity has only increased 
a factor of 1.7 over the initial value — clearly not enough 
to stabilize the equilibrium. Similar behavior is found 
for all of our runs that reach a cooling catastrophe. 

5.3. Dependence on Cluster Parameters 

We now consider the effect of varying the cluster pa- 
rameters, as exemplified by runs T2 and T3 in Table 
We will only compare the tangled magnetic field geome- 
tries since that is the most physically relevant set-up. 
Run T2 is modeled on Abell 2390, a hot, massive clus- 
ter. This NFW halo has a larger mass and scale radius. 
The core temperature rises from a central temperature 
of 4 keV to an outer temperature of 9.5 keV at 300 kpc. 
We choose the softening radius to be approximately 1/20 
of the scale radius as we did for our fiducial case. The 
dependence of our results on cluster mass and tempera- 
ture is small. The HBI growth time and central cooling 
time for our model of A2390 are 130 Myr and 1.1 Gyr, 
respectively, similar to our fiducial run. This run evolves 
in a similar way to run Tl, reaching the cooling catastro- 
phe in 3.2 Gyr with similar saturated parameters. Due 
to the vigor of the cooling catastrophe in this more mas- 
sive case, the run does not reach 9.5 Gyr, as most of the 
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Fig. 11. — The color scale shows the temperature in keV and the arrows represent the magnetic field unit vector in the x—y plane for 
run Rl with an initially radial field. The plots are at t = 0, 1.6, 4.8, and 9.5 Gyr (from left to right and top to bottom). As the magnetic 
field becomes more azimuthally wrapped, the cluster core reaches the cooling floor. 
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Fig. 12. — The shell-averaged magnetic field angle versus radius 
for run Tl-HB. At 2.2 Gyr, radial inflow has partially reoriented 
the field line radially at ~ 20 kpc. 

mass has collapsed to the central few zones. 
In run T3, we examine the effect of changing the ini- 



tial thermal profile but not the NFW parameters. This 
run is the same as the fiducial case but the temperature 
now initially varies from 1 to 6 keV. The physics of how 
the magnetic geometry is modified remains very similar; 
however, the cooling catastrophe occurs at an earlier time 
of 1.5 Gyr. The primary reason for this is that the higher 
initial central density and lower temperature make this 
cluster cool faster to the temperature floor. Thus, we see 
only a very weak dependence on the cluster parameters, 
mostly being driven by the initial location on the cooling 
curve. 

5.4. Dependence on Central Entropy 

Motivated by the discussion of the role of central en- 
tropy in we have undertaken a parameter study in 
central entropy. The runs Tl, El, E2, and E3 form a 
monotonic progression of central entropy from 22.4 to 
122 keV cm 2 . We have effected the entropy variation 
by modifying the initial central temperature while main- 
taining thermal equilibrium. We fit our cluster entropy 
profiles with a power-law of the form of equation Q. The 
parameter Ko is typically within 10-20% of the central 
entropy. The primary difference evident in examining 
Table [3] is that the time of the cooling catastrophe in- 
creases as Kq increases-reaching a maximum of just over 
10 Gyr for the highest entropy case with nG magnetic 
fields (run E3). This phenomenon is easy to understand 
since the central cooling time itself increases as Kq in- 
creases. 
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Fig. 13. — Azimuthally-averaged temperature profiles for our 
high entropy (85 keV cm 2 ) and high magnetic field (1 fiG) run (E2- 
HB). The HBl-driven cooling catastrophe is significantly stabilized. 
The cluster does not hit the cooling floor until 7.0 Gyr. 

It is interesting to compare these runs to several high- 
entropy simulations without conduction, runs E2-NC 
and E3-NC. These runs have cooling but no thermal con- 
duction. First, it is clear that the central cooling time 
estimated by equation (|2ip is an overestimate compared 
to the actual time of the cooling catastrophe. For exam- 
ple, run E3 has a predicted cooling timescale of 9.3 Gyr 
but in fact reaches a cooling catastrophe without conduc- 
tion in 5.9 Gyr. Second, these runs without conduction 
show that even though the effective conductivity is re- 
duced by the HBI, the time to a cooling catastrophe is 
significantly longer than in the absence of conduction. In 
the case of run E3, thermal conduction delays the cooling 
catastrophe from 5.3 Gyr to >9.5 Gyr. 

We now consider runs Tl-HB and E2-HB which have 
higher magnetic fields of 1 fxG. In the former, the low en- 
tropy run, magnetic tension does very little to stabilize 
the HBI resulting in a cooling catastrophe at almost the 
same time as the 1 nG case. In the latter, the high en- 
tropy run, magnetic tension plays a more significant role 
and increases the time of the cooling catastrophe from 
5.7 Gyr (the low field case) to 7.0 Gyr (the higher field 
case). Figure [T31 shows the evolution of the temperature 
profile for run E2-HB. The cooling is especially slow when 
the temperature is high. In fact, the time to reach the 
cooling catastrophe, ~ 7 Gyr, is longer than the typical 
time between major merg ers for galaxy clusters ~ 4-5 
Gyr dCohn fc Whitdl200l . 

How does a cluster reach the high entropy states for 
which the cooling catastrophe can be avoided? Cluster 
mergers may provide shocks that heat the cluster and 
boost its entropy. In addition, strong AGN feedback may 
increase the cluster entropy enough to slow the cooling 
instability, even in the presence of the HBI. 

5.5. Isothermal Initial Conditions 



As a further exploration of the important physics in 
clusters, we present the results of two simulations, runs 
Isol and Iso2, that are initialized with isothermal tem- 
perature profiles at 4.5 keV. Our usual approach of im- 
posing thermal equilibrium (eq. [23]) does not work in 
this case, and the central density becomes a free parame- 
ter for constructing the hydrostatic equilibrium. By def- 
inition, the initial state is not in thermal equilibrium. 
The other initial parameters are the same as our fiducial 
case. 

Run Isol has an initial central density of n e n = 2.5 x 
10~ 2 , which corresponds to an entropy of 52.5 keV cm' 1 
and a cooling time of 2.6 Gyr. By virtue of this short 
cooling time and the action of the HBI, the core expe- 
riences the cooling catastrophe at 2.1 Gyr. By lowering 
the central density by a factor of 5 to n e n = 5 x 1CP 3 , 
run Iso2 has a central entropy of ~ 154 keV cm 2 and a 
cooling time of 12.8 Gyr. After 9.5 Gyr, this run has 
developed a slightly relaxed core (Tj w 4.1 keV), but is 
very far from the cooling catastrophe. 

Starting from the non-equilibrium initial condition, we 
see very similar qualitative behavior to our equilibrium 
model. Runs with short cooling times develop both the 
cool core and the HBI quickly, and runs with long cooling 
times are much more stable. In our isothermal clusters 
the state is merely transient, as even clusters with long 
central cooling times and high entropy would eventually 
evolve into relaxed, cool-core clusters. 

5.6. Resolution Dependence 

We perform one high resolution simulation of our fidu- 
cial case, Tl-256, which is a (256) 3 tangled field simula- 
tion. We terminate this run at 6.3 Gyr (2/3 of the normal 
run time) in light of the large processor time required. 
All of the qualitative results discussed thus far hold up at 
higher resolution. We do see some minor differences re- 
sulting from the increased resolution: smaller scales are 
now available on which the HBI can act. The HBI both 
amplifies the magnetic energy slightly more and is able 
to reach an even larger average magnetic field angle. The 
final volume-averaged magnetic field angle of 77° corre- 
sponds to an incredibly azimuthally wrapped magnetic 
field with a tiny effective conductivity of /s p ~ 0.034. 
This precipitates the cooling catastrophe on a slightly 
shorter timescale. 

5.7. Experiments with Central Heating 

Surveys of galaxy cluster cores consistently find X- 
ray cavities or bubbl es filled with radio emitting plasma 
or co smic rays (e.g., iBlrzan et al]|2004t |D unn fc Fabian! 
2006). These structures indicative of feedback are es- 
pecially prevalent in clusters with short central cooling 
times, roughly the same population of low entropy clus- 
ters mentioned previously. In an effort to understand the 
interplay between the HBI and heating, we proceed with 
a preliminary analysis of heating in cluster cores. 

A number of groups have proposed cosmic rays from 
a central AGN as a heating mechanism. In partic- 
ular, streaming cosmic rays can excite Alfven waves 
which nonlinearly Landa u d amp to hea t the plasma 
(jLoewenstein et all 11991ft . IGuo fc Oh! (|2008ft have 
demonstrated in ID models that a combination of pa- 
rameterized cosmic ray feedback and conduction can pre- 
vent significant cooling for a Hubble time. For our test 
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Fig. 14. — The thermal evolution of our fiducial atmosphere with 
heating in the central 20 kpc, run HI. The cooling catastrophe is 
pushed outwards in radius. 

problems, our heating luminosity is parameterized as in 
equation (jTTJ) with the initial normalizat ion set by equa- 
tion (HHJ), motivated bv lChandranl (|2005f ). These heating 
functions are generic and do not discriminate among cos- 
mic ray or mechanical energy injection. 

The feedback dynamics of the cluster core is qualita- 
tively simple. As the core cools, the accretion rate onto 
the supermassive (~ 10 9 M Q ) black hole at the center of 
the cD galaxy increases. As M increases, the feedback 
heating increases, slowing the cooling. If heating be- 
comes too effective, the accretion ceases, and a feedback 
loop is established. Thus, a simple static heating model 
is insufficient, and we instead implement a rough time- 
variable version. Namely, we sum the cooling luminosity 
within the heating effective radius, r# , at t = 0, 



£n = 



47rr 3 /3 ' 



(33) 



We then calculate the cooling luminosity in a similar way 
at every timestep and scale the initial heating luminosity 
to the current cooling luminosity as 



H(t) 



£(t) 



C(t = 0) 



Ho. 



(34) 



This methodology is not ideal, but given that we cannot 
resolve the Bondi radius on our grid, it is preferable to ex- 
trapolating the central density and temperature down to 
the Bondi radius; and it ensures we have an approximate 
feedback mechanism. It should be noted that this heat- 
ing model can be numerically unstable when the cooling 
instability has progressed, and large feedback heating is 
added within a small region. We will improve this treat- 
ment in future work. As an example calculation, we take 
our fiducial atmosphere and add a total initial heating 
luminosity of Ltherm = 10 43 erg s _1 with a characteristic 
radius of r# = 20 kpc. We construct an initial thermal 



equilibrium such that heating, cooling, and thermal con- 
duction all balance. This run, labeled HI, has a very 
interesting thermal evolution as is shown in Figure [T4l 

The HBI proceeds very slowly for this simulation since 
the initial average magnetic field is 3.5/iG, strong enough 
to exert significant tension. What is especially interest- 
ing about this run is that the centrally concentrated feed- 
back heating drives a minimum in the temperature pro- 
file at ^20 kpc after 5.4 Gyr of evolution. This type of 
profile, with a minimum slightly offset from the center, 
is ac tually observed in some clusters (|Sanderson et al.1 
2006). The cooling flow region seems to be simply pushed 
outwards from the center of the cluster, perhaps evi- 
dence that an additional volumetric heating component 
is needed. The heating power at 6 Gyr has increased 
only modestly to a new value of 6.6 x 10 43 erg s -1 . Due 
to the combination of heating and the HBI the cooling 
catastrophe occurs much later than the initial central 
cooling time of 1.3 Gyr. Unfortunately, we are not able 
to follow this run to completion as the sharp tempera- 
ture discontinuity combined with the rapid cooling that 
follows leads to numerical instabilities. While not neces- 
sarily thermally stable, this run demonstrates that clus- 
ter cores with heating and the HBI can remain stable for 
longer than the time between major mergers, although 
not necessarily a Hubble time. In general, it appears 
that the combination of magnetic fields of > 1 /iG (to 
slow the HBI) and a modest amount of AGN feedback 
significantly slow the cooling catastrophe, even for low- 
entropy clusters. 

6. DISCUSSION AND CONCLUSIONS 

The plasma in the intracluster medium of galaxy clus- 
ters is dilute and magnetized with a mean free path large 
compared to the gyroradius. Under these conditions, 
heat transport is anisotropic along magnetic fields. This 
results in the ICM being unstable to the heat-flux-driven 
buoyancy instabili ty in regions w here the temperature 
increases outward (|Quataertll2008( ). The cores of galaxy 
clusters also often have short cooling times of < 500 Myr. 
In order to understand the thermal evolution of galaxy 
clusters with cooling and the HBI, we have performed 
three-dimensional time-dependent MHD simulations of 
galaxy clusters cores. 

Isolated galaxy clusters evolved with magnetic fields, 
anisotropic conduction, and cooling share a number of 
common properties. We begin with a cluster that is in 
both hydrostatic and thermal equilibrium. After ~ 100 
Myr, the HBI begins to rearrange the magnetic field ge- 
ometry in the cluster core; the magnetic field saturates 
with an average angle between the magnetic field and the 
radial direction of ~75°. Second, as the magnetic geom- 
etry is rearranged to be tangential to the temperature 
gradient, the magnetic field exerts a thermally insulating 
effect, reducing the effective radial thermal conductivity 
to < 10% of the Spitzer value. Finally, having reduced 
the thermal conduction from the outer parts of the clus- 
ter, and lacking another heat source, the core proceeds 
to a cooling catastrophe on a timescale comparable to 
the initial central cooling time. 

We have studied a number of different parameter vari- 
ations and find several in teresting trends. Motivated by 
the observational work of lVoit et al.l (|2008f ). we have ex- 
plored the effects of different initial entropies. For larger 
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central cluster entropies, the time of the cooling catas- 
trophe is delayed to more than 9.5 Gyr for our highest 
entropy cluster (122 keV cm 2 ). In addition, we find that 
stronger magnetic fields, > 3.5 /iG, can suppress the HBI 
via magnetic tension forces. The onset of the cooling 
catastrophe can thus be delayed in these higher magnetic 
field calculations. 

We have also carried out initial calculations of the ef- 
fects of heating on the ICM using a parameterized heat- 
ing function in which the total heating power is propor- 
tional to the total rate of cooling in the central 20 kpc. 
Despite some difficulty with numerical instabilities inher- 
ent in the method, we find that modest heating rates of 
10 43 erg s _1 can substantially delay the cooling catas- 
trophe to > 5 Gyr, longer than the time between ma- 
jor cluster mergers. There is, however, some evidence 
that centrally concentrated heating may simply move the 
cooling catastrophe further out in the core (Fig. [T"4"]) . 

Observations and simplified theoretical models both 
suggest that t here are t wo d is tinct quas i- stable 
cluster states dVoit et all 120081: iGuo fc Ohl 120081 : 
ICavagnolo et al.l l2009f ) . High entropy, fairly isothermal 
cluster cores have long growth times for both the Field 
instability and the HBI. The longer central cooling times 
require less conductive heating to balance cooling. These 
thermal states are long-lived even in the absence of merg- 
ers. By contrast, low entropy relaxed clusters (cool-core 
clusters) with short central cooling times are unstable 
to both the Field instability and the HBI. As we have 
shown, this population of clusters cannot be stabilized 
by conduction alone and must have an additional feed- 
back mechanism, plausibly the central AGN but poten- 
tially other sources. Observations of cool core clusters 
with central entropies Kq < 30 keV cm 2 show a number 
of feedback indicators, including Ha emission indicative 
of cool gas at 10 4 K, radio emission indicative of AGN 
feedback, and optical color grad ients indicative of cen- 
tral star forma tion in the BCG (|Cavagnolo et al.l 120081 : 
IVoit et all [2001 ). High entropy clusters, in which con- 
duction is more important, generally show none of these 
feedback indicators. Our simulations of these high en- 
tropy clusters show that they are thermally stable for 
cosmologically long timescales, and that conduction pro- 
vides a significant stabilizing effect, e.g. runs E2 and E3 
(see, Tables 2 and 3). 

It may be possible for clusters to transition between 
these two populations. Relaxed clusters may be pro- 
moted to high entropy clusters by significant heating, 
such as a major merger or an especially energetic feed- 
back event. Recent work shows that disruption of cool 
cores in a merger is possible at cosmologically early 
times but difficult at late times (|Burns "eFld1l2008h . Al- 
ternatively, isothermal moderate entropy clusters can 
eventually become relaxed cool-core clusters over long 
timescales. 
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A key task for future work is to better understand the 
proposed heating mechanisms for low entropy clusters. 
In particular, bubbles and jets from an AGN are far 
from geometrically isotropic. Not only must the heat- 
ing be locally efficient, but the heating must then be 
distributed by some mechanism azimuthally around the 
cluster core to prevent a cooling catastrophe. The en- 
hanced azimuthal heat transport from the HBI may play 
a significant role in redistributing local AGN heating 
throughout the cluster core. 

In future work, we will examine these heating mecha- 
nisms in more detail including the relevant physics. For 
the case of buoyant bubbles, there are many unanswered 
questions about the disruption time of the bubbles. This 
shredding is governed by Rayleigh- Taylor and Kelvin- 
Helmholtz instabilities. In the full Braginskii-MHD 
treatment, momentum is transported by ions anisotropi- 
cally along magnetic field lines. If the bubbles are indeed 
draped by magnetic fields, then the RT and KH instabil- 
ities will be modified by an anisotropic Braginskii viscos- 
ity. Cosmic rays may play a role in directly heatin g the 
plasma by exciting Alfven waves (|Guo fc Obl l2008). Fi- 
nally, galaxy wakes in a full co smological context can als o 
provide heating to the ICM (|Conrov fc Ostrikerl [20081) . 
They may also increase the importance of thermal con- 
duction by competing with the HBI to reorient the mag- 
netic field. 

A key lesson of this work is that it is difficult to char- 
acterize the ICM plasma as having a single thermal con- 
ductivity parameterized by a constant /s p . Buoyancy 
instabilities such as the HBI and MTI directly modify 
the magnetic geometry and self-consistently evolve the 
system to a new state that may enhance or suppress the 
effective conductivity. In the cores of galaxy clusters, the 
HBI suppresses thermal conduction from the large heat 
reservoir at large radii. In the absence of AGN feedback 
or very large magnetic fields in cluster cores, it appears 
that conduction alone cannot solve the cooling flow prob- 
lem. 
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